########################################################################################################################################
## ECONOMIC SHOCKS & ISSUE IMPORTANCE                                                                                                 ##
## Author:  Sarah Gomm, based on Hanretty et al. 2020 https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/0EIQXL ##
## Date:    February 2023                                                                                                             ##
## Data:    SEP Wave 6                                                                                                                ##
########################################################################################################################################


rm(list=ls())
# Check the username and set the working directory   ---------------------------
# Get the username 
username <- Sys.getenv("USERNAME")

if (username == "FILL IN YOUR USERNAME") {
  setwd("FILL IN YOUR WORKING DIRECTORY")
} else {
  stop("Unknown user. Please check your username.")
}



#Install & load packages   -----------------------------------------------------
packages <- c("ggplot2", "knitr", "pander", "reshape2", "xtable") 

for (p in packages) {
  if (p %in% installed.packages()[,1]) require(p, character.only=T)
  else {
    install.packages(p)
    library(p, character.only=T)
  }
  if (p %in% rownames(.packages()) == FALSE)
  {library(p, character.only=T)
  }
}


# ------------------------------------------------------------------------------
#  Issue importance for full sample  (Appendix)     ----------------------------
# ------------------------------------------------------------------------------

#Load data     -----------------------------------------------------------------
load("Surveydata_fullsample.Rdata")
load("Posterior_fullsample.Rdata")


#Data preperation: full sample  (cf. Hanretty et al. 2020)    ------------------
N_items <- 16
N_respondents <- nrow(svyclean)
differing_positions_count <- apply(conjoint_scaling_data$candidatepositionID,c(1,4),function(x) sum(x[1,] != x[2,]))

# cutpoints for A / B
gamma_est <- apply(rstan::extract(posterior)$gamma,c(2),mean)
print(round(gamma_est,2))

# positions of alternatives
psi_est <- apply(rstan::extract(posterior)$psi,c(2,3),mean)
psi_int <- apply(rstan::extract(posterior)$psi,c(2,3),quantile,c(0.025,0.975))
rownames(psi_est) <- issue_labels
colnames(psi_est) <- c("Disagree", "Rather Disagree", "Rather Agree", "Agree")
print(round(psi_est,2))

psi_sd <- apply(rstan::extract(posterior)$psi,c(2,3),sd)

# check ordering
order_check <- apply(psi_est,1,diff)

# check posterior probability of correct ordering
psi_diff_chain <- apply(rstan::extract(posterior)$psi,c(1,2),diff)
psi_diff_p <- apply(psi_diff_chain,c(1,3),function(x) mean(x < 0))

# population proportions for each position
pi_est <- apply(rstan::extract(posterior)$pi,c(2,3),mean)
pi_se <- apply(rstan::extract(posterior)$pi,c(2,3),sd)
rownames(pi_est) <- issue_labels

chi_chain <- rstan::extract(posterior)$chi
chi_est <- apply(chi_chain,c(2),mean)
chi_int <- t(apply(chi_chain,c(2),quantile,c(0.025,0.975)))
names(chi_est) <- rownames(chi_int) <- issue_labels
round(chi_int,2)

chi_compare_chain <- array(NA,c(nrow(chi_chain),N_items,N_items))
for (i in 1:nrow(chi_chain)) chi_compare_chain[i,,] <- outer(chi_chain[i,],chi_chain[i,],FUN=">")
chi_compare_p <- apply(chi_compare_chain,c(2,3),mean)
rownames(chi_compare_p) <- colnames(chi_compare_p) <- issue_labels

disagreement_table <- data.frame(issue=issue_labels,
                                 disagreement=chi_est,
                                 disagreement_lb=chi_int[,1],disagreement_ub=chi_int[,2])
sort_order <- order(disagreement_table$disagreement,decreasing=TRUE)
inv_sort_order <- order(disagreement_table$disagreement,decreasing=FALSE)
disagreement_league_table <- disagreement_table[sort_order,]
rownames(disagreement_league_table) <- NULL

support_table <- pi_est[sort_order,]

responses_by_differing_count <- prop.table(table(cbind(as.numeric(differing_positions_count)),as.numeric(conjoint_scaling_data$Y)),1)

#Alternative Support Estimates
colnames(support_table) <- c("Disagree", "Rather Disagree", "Rather Agree", "Agree")

print(xtable(support_table,
             caption = "Population weighted support for each of the four policy alternatives on each issue, based on self-reported response questions.",
             digits=2
),
type="latex",
table.placement="!h", comment=FALSE, file = "Appendix_Table_A12.tex", include.rownames = FALSE
)

#Alternative Position estimates
round2 <- function(x) format(round(x,2),nsmall=2)
position_table <- round2(psi_est)
colnames(position_table) <- c("Disagree", "Rather Disagree", "Rather Agree", "Agree")

position_table_se <- apply(round2(psi_sd),2,function(x) paste0("(",x,")"))
colnames(position_table_se) <- rep("SE",4)

position_table <- cbind(position_table,position_table_se)
position_table <- position_table[sort_order,c(1,5,2,6,3,7,4,8)]

print(xtable(position_table,
             caption = "Posterior mean relative positions of the four policy alternatives on each issue. Posterior standard deviations (SE) in parentheses.",
             digits=2
),
type="latex",
table.placement="!h", comment=FALSE, sanitize.colnames.function=NULL, file = "Appendix_Table_A13.tex"
)

# Table with Importance Estimates for full sample ------------------------------
disagreement_league_table_print <- data.frame(
  Importance = round2(disagreement_league_table$disagreement),
  Interval = paste0("(",round2(disagreement_league_table$disagreement_lb),"-",round2(disagreement_league_table$disagreement_ub),")")
)
rownames(disagreement_league_table_print) <- disagreement_league_table$issue

print(xtable(disagreement_league_table_print,
             caption = "Posterior mean importance for each issue, with central 95\\% posterior interval in parentheses.",
             digits=2
),
type="latex",
table.placement="!h", comment=FALSE, sanitize.colnames.function=NULL, file = "Appendix_Table_A14.tex"
)


# ------------------------------------------------------------------------------
# Issue importance for control group only  (Main Article)  ---------------------
# ------------------------------------------------------------------------------

rm(list= ls()[!(ls() %in% c('round2','issue_labels','issue_labels_signed', 'N_items'))])
load("Surveydata_control.Rdata")
load("Posterior_control.Rdata")

#Data preparation for control group only
differing_positions_count <- apply(conjoint_scaling_data$candidatepositionID,c(1,4),function(x) sum(x[1,] != x[2,]))

# cutpoints for A / B
gamma_est <- apply(rstan::extract(posterior)$gamma,c(2),mean)
print(round(gamma_est,2))

# positions of alternatives
psi_est <- apply(rstan::extract(posterior)$psi,c(2,3),mean)#Betrag von!
psi_int <- apply(rstan::extract(posterior)$psi,c(2,3),quantile,c(0.025,0.975))
rownames(psi_est) <- issue_labels
colnames(psi_est) <- c("Disagree", "Rather Disagree", "Rather Agree", "Agree")
print(round(psi_est,2))

psi_sd <- apply(rstan::extract(posterior)$psi,c(2,3),sd)

# check ordering
order_check <- apply(psi_est,1,diff)

# check posterior probability of correct ordering
psi_diff_chain <- apply(rstan::extract(posterior)$psi,c(1,2),diff)
psi_diff_p <- apply(psi_diff_chain,c(1,3),function(x) mean(x < 0))

# population proportions for each position
pi_est <- apply(rstan::extract(posterior)$pi,c(2,3),mean)
pi_se <- apply(rstan::extract(posterior)$pi,c(2,3),sd)
rownames(pi_est) <- issue_labels

chi_chain <- rstan::extract(posterior)$chi
chi_est <- apply(chi_chain,c(2),mean)
chi_int <- t(apply(chi_chain,c(2),quantile,c(0.025,0.975)))
names(chi_est) <- rownames(chi_int) <- issue_labels
round(chi_int,2)

chi_compare_chain <- array(NA,c(nrow(chi_chain),N_items,N_items))
for (i in 1:nrow(chi_chain)) chi_compare_chain[i,,] <- outer(chi_chain[i,],chi_chain[i,],FUN=">")
chi_compare_p <- apply(chi_compare_chain,c(2,3),mean)
rownames(chi_compare_p) <- colnames(chi_compare_p) <- issue_labels

disagreement_table <- data.frame(issue=issue_labels,
                                 disagreement=chi_est,
                                 disagreement_lb=chi_int[,1],disagreement_ub=chi_int[,2])
sort_order <- order(disagreement_table$disagreement,decreasing=TRUE)
inv_sort_order <- order(disagreement_table$disagreement,decreasing=FALSE)
disagreement_league_table <- disagreement_table[sort_order,]
rownames(disagreement_league_table) <- NULL

support_table <- pi_est[sort_order,]

responses_by_differing_count <- prop.table(table(cbind(as.numeric(differing_positions_count)),as.numeric(conjoint_scaling_data$Y)),1)

#Alternative Support Estimates for control group
colnames(support_table) <- c("Disagree", "Rather Disagree", "Rather Agree", "Agree")
support_table 

#Alternative Position estimates for control group 
round2 <- function(x) format(round(x,2),nsmall=2)
position_table <- round2(psi_est)
colnames(position_table) <- c("Disagree", "Rather Disagree", "Rather Agree", "Agree")

position_table_se <- apply(round2(psi_sd),2,function(x) paste0("(",x,")"))
colnames(position_table_se) <- rep("SE",4)

position_table <- cbind(position_table,position_table_se)
position_table <- position_table[sort_order,c(1,5,2,6,3,7,4,8)]
position_table

# Table with Importance Estimates for control group ----------------------------
disagreement_league_table_print <- data.frame(
  Importance = round2(disagreement_league_table$disagreement),
  Interval = paste0("(",round2(disagreement_league_table$disagreement_lb),"-",round2(disagreement_league_table$disagreement_ub),")")
)
rownames(disagreement_league_table_print) <- disagreement_league_table$issue
disagreement_league_table_print

# ------------------------------------------------------------------------------
# Figure 1 (main article): Issue importance for control group 
# ------------------------------------------------------------------------------

#Main Plot for issue importance in control group 
labs = issue_labels[inv_sort_order]
sublabs_green = c("Moratorium gen. plants", "Animal-friendly care", "Protection of predators", "Construction of buildings")
sublabs_econ = c("Financial support childcare", "Minimum wage of 4000CHF", "Opening hours of retailers", "Raising retirement age to 67")
sublabs_econgreen = c("Renounce fossil use 2050", "Extenstion of CO2 charge", "Six lanes for motorways", "Subsidies for all farmers")


tiff("Main_Figure1.tiff", units="in", width=9, height=6, res=300)

m <- rbind(c(1,1,1,2,2))
layout(m)

par(mar=c(5,14,5,2))
plot(0,0,type="n",xlim=c(0,3.5),ylim=c(0,N_items+1),axes=FALSE,xlab="Estimated Policy Alternative Locations",ylab="", cex.lab = 1.4)
axis(2,las=2,at=1:16,labels=issue_labels[inv_sort_order],cex.axis=1.1) 
axis(2,at=match(sublabs_green,labs),labels=sublabs_green,las=2, adj=1, cex.axis=1.1, col.axis="blue")
axis(2,at=match(sublabs_econ,labs),labels=sublabs_econ,las=2, adj=1, cex.axis=1.1, col.axis="darkred")
axis(2,at=match(sublabs_econgreen,labs),labels=sublabs_econgreen,las=2, adj=1, cex.axis=1.1, col.axis="lightblue")
axis(1)
for (j in 1:N_items){
  points(psi_est[inv_sort_order[j],],rep(j,4),pch=16,cex=6.2*sqrt(pi_est[inv_sort_order[j],]),col=rgb(0,0,0,0.2))
  text(psi_est[inv_sort_order[j],],rep(j,4),1:4,cex=3.2*sqrt(pi_est[inv_sort_order[j],]),col=rgb(0,0,0,0.5))
  text(3.5,j,round2(disagreement_league_table[N_items-j+1,2]),cex=1.1)
}

par(mar=c(5,0,5,14))
plot(0,0,type="n",xlim=c(0,1.5),ylim=c(0,N_items+1),axes=FALSE,xlab="Estimated Importance",ylab="", cex.lab = 1.4)
axis(4,las=2,at=1:16,labels=issue_labels[inv_sort_order],cex.axis=1.1)
axis(4,at=match(sublabs_green,labs),labels=sublabs_green,las=2, adj=1, cex.axis=1.1, col.axis="blue")
axis(4,at=match(sublabs_econ,labs),labels=sublabs_econ,las=2, adj=1, cex.axis=1.1, col.axis="darkred")
axis(4,at=match(sublabs_econgreen,labs),labels=sublabs_econgreen,las=2, adj=1, cex.axis=1.1, col.axis="lightblue")
axis(1)
for (j in 1:N_items){
  lines(disagreement_league_table[N_items-j+1,3:4],rep(j-0.25,2),col=rgb(0,0,0,0.5),lwd=2)
  text(disagreement_league_table[N_items-j+1,2],j,round2(disagreement_league_table[N_items-j+1,2]),cex=1.1)
}

dev.off()



# ------------------------------------------------------------------------------
# Issue importance as a function of economic affectedness (Main Article) -------
# ------------------------------------------------------------------------------

rm(list= ls()[!(ls() %in% c('round2','issue_labels','issue_labels_signed', 'N_items', 'inv_sort_order'))])

# Table for issue importance dependent on  worse work situation ----------------
load("Surveydata_fullsample.Rdata")
load("Posterior_fullsample_by-treatment.Rdata")

beta_chain <- rstan::extract(posterior)$beta 
N_sim <- dim(beta_chain)[1]
beta_est <- apply(beta_chain,c(2,3),mean)
beta_se <- apply(beta_chain,c(2,3),sd)
rownames(beta_est) <- issue_labels
colnames(beta_est) <- c( "Treatment (economic shock)")

# calculate which coefficients are significantly different from the mean coefficient for that variable
beta_mean_chain <- apply(beta_chain,c(1,3),mean)
beta_sign_chain <- array(NA,dim=dim(beta_chain))
for (sim in 1:N_sim) for (j in 1:N_items) beta_sign_chain[sim,j,] <- beta_chain[sim,j,] > beta_mean_chain[sim,]
beta_p <- apply(beta_sign_chain,c(2,3),mean)
beta_sig <- -1*(beta_p < 0.05) + (beta_p > 0.95)
beta_star <- array(c("(-)","","(+)")[beta_sig+2],dim(beta_sig))
beta_print <- data.frame(round(beta_est,3),beta_star)
beta_print <- beta_print[rev(inv_sort_order),]

#Confidence intervals beta coefficient
beta_est <- apply(beta_chain,c(2,3),mean)
beta_int <- t(apply(beta_chain,c(2),quantile,c(0.025,0.975)))

beta_treatment_teable <- data.frame(issue=issue_labels,
                                    beta_treatment=beta_est,
                                    beta_treatment_lb=beta_int[,1],beta_treatment_ub=beta_int[,2])

beta_treatment_table_print <- data.frame(
  Beta_Treatment = round2(beta_treatment_teable$beta_treatment),
  Interval = paste0("(",round2(beta_treatment_teable$beta_treatment_lb),"-",round2(beta_treatment_teable$beta_treatment_ub),")"),
  Significance = beta_star
)
rownames(beta_treatment_table_print) <- beta_treatment_teable$issue

# add average across issues
beta_treatment_table_print[N_items+1,c(2)] <- ""
beta_treatment_table_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_treatment_table_print)[N_items+1] <- "Average"
colnames(beta_treatment_table_print)[c(3)] <- ""


# For interpretation: Calculate multiplicative effect from importance ----------
beta_multi <- -1*(1-exp(beta_est))
percent_matrix <- matrix(NA, nrow = nrow(beta_multi), ncol = ncol(beta_multi))

for (i in 1:nrow(beta_multi)) {
  for (j in 1:ncol(beta_multi)) {
    percent_matrix[i, j] <- paste0(signif(beta_multi[i, j] * 100, digits = 2), "%")
  }
}

multi_average <- paste0(round((-1 * (1 - exp(-0.107))) * 100, digits = 1), "%")
percent_matrix <- c(percent_matrix, multi_average)

# add the percent_matrix as a new column to the data frame
df <- cbind(beta_treatment_table_print, percent_matrix)
colnames(df) <- c( "Beta coefficients for treatment", "Interval", "", "Multiplicative effect")

# Print Table for Appendix (A.15)
print(xtable(df,
             caption = "Coefficient estimates for variation in importance given a worsened work situation. Coefficients
significantly higher than or lower than the average coefficient for the treatment variable across all issues are marked
with (+) or (-), respectively. The multiplicative effect is the natural exponent of the beta coefficient, indicating the impact of being treated on the disutility weight attached to voter-candidate disagreements."
             ,digits=2),
      type="latex",
      table.placement="!h", comment=FALSE, file = "Appendix_Table_A15.tex"
)


# Plot for issue importance dependent on economic affectedness -----------------

# Cutpoints for A / B
gamma_est <- apply(rstan::extract(posterior)$gamma,c(2),mean)

# exponent for utility loss function
beta_est <- apply(rstan::extract(posterior)$beta,c(2,3),mean)
beta_se <- apply(rstan::extract(posterior)$beta,c(2,3),sd)
rownames(beta_est) <- issue_labels
colnames(beta_est) <- colnames(conjoint_scaling_data$X)

# positions of alternatives
psi_est <- apply(rstan::extract(posterior)$psi,c(2,3),mean)
rownames(psi_est) <- issue_labels
colnames(psi_est) <- c("Disagree", "Rather Disagree", "Rather Agree", "Agree")

# population proportions for each position
pi_est <- apply(rstan::extract(posterior)$pi,c(2,3),mean)
rownames(pi_est) <- issue_labels

# calculate importance by treatment
ordinal_response_table <- matrix(NA,N_items,4)
for (j in 1:N_items) ordinal_response_table[j,] <- prop.table(xtabs(svyclean$weight~ordinal_responses[,j]))


# Confidence intervals issue importance   --------------------------------------

# Calculate Confidence intervals for beta 
confidence_level <- 0.95
beta_est_lower <- apply(beta_chain, c(2, 3), function(x) quantile(x, (1 - confidence_level) / 2))
beta_est_upper <- apply(beta_chain, c(2, 3), function(x) quantile(x, 1 - (1 - confidence_level) / 2))

## positions of alternatives
psi_est_lower <- apply(rstan::extract(posterior)$psi, c(2, 3), function(x) quantile(x, 0.025))
psi_est_upper <- apply(rstan::extract(posterior)$psi, c(2, 3), function(x) quantile(x, 0.975))


# population proportions for each position
pi_est <- apply(rstan::extract(posterior)$pi,c(2,3),mean)
rownames(pi_est) <- issue_labels

# calculate importance by worse work situation for lower bound
ordinal_response_table <- matrix(NA,N_items,4)
for (j in 1:N_items) ordinal_response_table[j,] <- prop.table(xtabs(svyclean$weight~ordinal_responses[,j]))

# calculate total weight of importance for each issue, by level of worse work situation
importance_lower <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est_lower[j,] * (g-1)))
    for (p in 1:4) importance_lower[j,g] <- importance_lower[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est_lower[j,p]-psi_est_lower[j,])))
  }
}


# calculate importance by worse work situation for upper bound
ordinal_response_table <- matrix(NA,N_items,4)
for (j in 1:N_items) ordinal_response_table[j,] <- prop.table(xtabs(svyclean$weight~ordinal_responses[,j]))

# calculate total weight of importance for each issue, by level of worse work situation
importance_upper <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est_upper[j,] * (g-1)))
    for (p in 1:4) importance_upper[j,g] <- importance_upper[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est_upper[j,p]-psi_est_upper[j,])))
  }
}


#Mean importance
importance <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est[j,] * (g-1)))
    for (p in 1:4) importance[j,g] <- importance[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est[j,p]-psi_est[j,])))
  }
}

rownames(importance) <- issue_labels

importance_table <- data.frame(issue=issue_labels,
                               importance_control=importance[,1],
                               importance_lower_control=importance_lower[,1],importance_upper_control=importance_upper[,1],
                               importance_treat=importance[,2],
                               importance_lower_treat=importance_lower[,2],importance_upper_treat=importance_upper[,2])


importance_table_print <- data.frame(
  Importance_Control = round2(importance_table$importance_control),
  Interval_control = paste0("(",round2(importance_table$importance_lower_control),"-",round2(importance_table$importance_upper_control),")"),
  Importance_Treat = round2(importance_table$importance_treat),
  Interval_Treat = paste0("(",round2(importance_table$importance_lower_treat),"-",round2(importance_table$importance_upper_treat),")")
)

rownames(importance_table_print) <- beta_treatment_teable$issue
colnames(importance_table_print) <- c("Control", " ", "Treatment", " ")

#app_imp_treat_betas.tex (in Latex file)
print(xtable(importance_table_print,
             caption = "Importance estimates given a worsened work situation or no change wit 95% posterior intervals."
             ,digits=3),
      type="latex",
      table.placement="!h", comment=FALSE, file = "Appendix_Table_A16.tex"
)


#Collapse table
imp_env <- data.frame(importance[c(1:4),])
imp_ecoenv <- data.frame(importance[c(5:8),])
imp_econ <- data.frame(importance[c(9:12),])
imp_oth <- data.frame(importance[c(13:16),])

# Calculate the means for each variable in each data frame
means_imp <- lapply(list(imp_env, imp_ecoenv, imp_econ, imp_oth), function(df) colMeans(df))

# Merge the means data frames into one data frame
merged_df <- do.call(rbind, means_imp)
rownames(merged_df) <- c("Environment", "Environment/Economy", "Economy", "Other")

#app_imp_treat_control_collapsed.tex (in Latex file)
print(xtable(merged_df,
             caption = "Importance estimates for treatment and control group."
             ,digits=3),
      type="latex",
      table.placement="!h", comment=FALSE, file = "Appendix_Table_A17.tex", include.rownames = FALSE
)

# ------------------------------------------------------------------------------
# Figure 2 (main article): Issue importance dependent on treatment -------------
# ------------------------------------------------------------------------------

tiff("Main_Figure2.tiff", units="in", width=22, height=14, res=600)

par(mfrow=c(2,2), mar=c(2, 0, 6, 0))
alpha_value <- 0.5

plot(0,0,xlim=c(-1,2),ylim=c(0,1.5),type="n",axes=FALSE,xlab="",ylab="", main="Environment-related issues", cex.main=2.8)
axis(1,at=c(0,1),labels=c("Not affected","Affected"),cex.axis=1.8)
for (j in 1:N_items) {
  lines(0:1,imp_env[j,],col="black", lwd = 2)
  lines(0:1,imp_ecoenv[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_econ[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_oth[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_env[j,1],rownames(imp_env)[j],pos=2,cex=1.8,col="black")
  text(1,imp_env[j,2],rownames(imp_env)[j],pos=4,cex=1.8,col="black")
  text(0,imp_ecoenv[j,1],rownames(imp_ecoenv)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_ecoenv[j,2],rownames(imp_ecoenv)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_econ[j,1],rownames(imp_econ)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_econ[j,2],rownames(imp_econ)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_oth[j,1],rownames(imp_oth)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_oth[j,2],rownames(imp_oth)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
}

plot(0,0,xlim=c(-1,2),ylim=c(0,1.5),type="n",axes=FALSE,xlab="",ylab="", main="Environment and economy-related issues", cex.main=2.8)
axis(1,at=c(0,1),labels=c("Not affected","Affected"),cex.axis=1.8)
for (j in 1:N_items) {
  lines(0:1,imp_env[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_ecoenv[j,],col="black", lwd = 2)
  lines(0:1,imp_econ[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_oth[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_env[j,1],rownames(imp_env)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_env[j,2],rownames(imp_env)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_ecoenv[j,1],rownames(imp_ecoenv)[j],pos=2,cex=1.8,col = "black")
  text(1,imp_ecoenv[j,2],rownames(imp_ecoenv)[j],pos=4,cex=1.8,col = "black")
  text(0,imp_econ[j,1],rownames(imp_econ)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_econ[j,2],rownames(imp_econ)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_oth[j,1],rownames(imp_oth)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_oth[j,2],rownames(imp_oth)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
}

plot(0,0,xlim=c(-1,2),ylim=c(0,1.5),type="n",axes=FALSE,xlab="",ylab="", main="Economy-related issues", cex.main=2.8)
axis(1,at=c(0,1),labels=c("Not affected","Affected"),cex.axis=1.8)
for (j in 1:N_items) {
  lines(0:1,imp_env[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_ecoenv[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_econ[j,],col="black", lwd = 2)
  lines(0:1,imp_oth[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_env[j,1],rownames(imp_env)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_env[j,2],rownames(imp_env)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_ecoenv[j,1],rownames(imp_ecoenv)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_ecoenv[j,2],rownames(imp_ecoenv)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_econ[j,1],rownames(imp_econ)[j],pos=2,cex=1.8,col="black")
  text(1,imp_econ[j,2],rownames(imp_econ)[j],pos=4,cex=1.8,col="black")
  text(0,imp_oth[j,1],rownames(imp_oth)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_oth[j,2],rownames(imp_oth)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
}


plot(0,0,xlim=c(-1,2),ylim=c(0,1.5),type="n",axes=FALSE,xlab="",ylab="", main="Other issues", cex.main=2.8)
axis(1,at=c(0,1),labels=c("Not affected","Affected"),cex.axis=1.8)
for (j in 1:N_items) {
  lines(0:1,imp_env[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_ecoenv[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_econ[j,],col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  lines(0:1,imp_oth[j,],col="black", lwd = 2)
  text(0,imp_env[j,1],rownames(imp_env)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_env[j,2],rownames(imp_env)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_ecoenv[j,1],rownames(imp_ecoenv)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_ecoenv[j,2],rownames(imp_ecoenv)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_econ[j,1],rownames(imp_econ)[j],pos=2,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(1,imp_econ[j,2],rownames(imp_econ)[j],pos=4,cex=1.8,col = rgb(0.8, 0.8, 0.8, alpha = alpha_value))
  text(0,imp_oth[j,1],rownames(imp_oth)[j],pos=2,cex=1.8,col="black")
  text(1,imp_oth[j,2],rownames(imp_oth)[j],pos=4,cex=1.8,col="black")
}

dev.off()



# ------------------------------------------------------------------------------
# SVP-voters: issue importance   -----------------------------------------------
# ------------------------------------------------------------------------------
#Issue Importance for SVP-voters
load("Surveydata_fullsample.Rdata")
load("Posterior_fullsample_by-partysvp.Rdata")

beta_chain <- rstan::extract(posterior)$beta 
N_sim <- dim(beta_chain)[1]
beta_est <- apply(beta_chain,c(2,3),mean)
beta_se <- apply(beta_chain,c(2,3),sd)
rownames(beta_est) <- issue_labels
colnames(beta_est) <- c( "Treatment (economic shock)")

# calculate which coefficients are significantly different from the mean coefficient for that variable
beta_mean_chain <- apply(beta_chain,c(1,3),mean)
beta_sign_chain <- array(NA,dim=dim(beta_chain))
for (sim in 1:N_sim) for (j in 1:N_items) beta_sign_chain[sim,j,] <- beta_chain[sim,j,] > beta_mean_chain[sim,]
beta_p <- apply(beta_sign_chain,c(2,3),mean)
beta_sig <- -1*(beta_p < 0.05) + (beta_p > 0.95)
beta_star <- array(c("(-)","","(+)")[beta_sig+2],dim(beta_sig))
beta_print <- data.frame(round(beta_est,3),beta_star)
beta_print <- beta_print[rev(inv_sort_order),]

#Confidence intervals beta coefficient
beta_est <- apply(beta_chain,c(2,3),mean)
beta_int <- t(apply(beta_chain,c(2),quantile,c(0.025,0.975)))

beta_treatment_teable <- data.frame(issue=issue_labels,
                                    beta_treatment=beta_est,
                                    beta_treatment_lb=beta_int[,1],beta_treatment_ub=beta_int[,2])

beta_treatment_table_print <- data.frame(
  Beta_Treatment = round2(beta_treatment_teable$beta_treatment),
  Interval = paste0("(",round2(beta_treatment_teable$beta_treatment_lb),"-",round2(beta_treatment_teable$beta_treatment_ub),")"),
  Significance = beta_star
)
rownames(beta_treatment_table_print) <- beta_treatment_teable$issue

# add average across issues
beta_treatment_table_print[N_items+1,c(2)] <- ""
beta_treatment_table_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_treatment_table_print)[N_items+1] <- "Average"
colnames(beta_treatment_table_print)[c(3)] <- ""

# add average across issues
beta_print[N_items+1,c(2)] <- ""
beta_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_print)[N_items+1] <- "Average"
colnames(beta_print)[c(2)] <- ""


# For interpretation: Calculate multiplicative effect from importance ----------
beta_multi <- -1*(1-exp(beta_est))
percent_matrix <- matrix(NA, nrow = nrow(beta_multi), ncol = ncol(beta_multi))

for (i in 1:nrow(beta_multi)) {
  for (j in 1:ncol(beta_multi)) {
    percent_matrix[i, j] <- paste0(signif(beta_multi[i, j] * 100, digits = 2), "%")
  }
}

multi_average <- paste0(round((-1 * (1 - exp(-0.077))) * 100, digits = 1), "%")
percent_matrix <- c(percent_matrix, multi_average)


# add the percent_matrix as a new column to the data frame
df <- cbind(beta_treatment_table_print, percent_matrix)
colnames(df) <- c( "Beta coefficients for treatment", "Interval", "", "Multiplicative effect")

#app_imp_treat_betas_multi.tex (in Latex file)
print(xtable(df,
             caption = ""
             ,digits=2),
      type="latex",
      table.placement="!h", comment=FALSE, file = "Appendix_Table_A18.tex"
)

# Calculate Confidence intervals for beta 
confidence_level <- 0.95
beta_est_lower <- apply(beta_chain, c(2, 3), function(x) quantile(x, (1 - confidence_level) / 2))
beta_est_upper <- apply(beta_chain, c(2, 3), function(x) quantile(x, 1 - (1 - confidence_level) / 2))

## positions of alternatives
psi_est_lower <- apply(rstan::extract(posterior)$psi, c(2, 3), function(x) quantile(x, 0.025))
psi_est_upper <- apply(rstan::extract(posterior)$psi, c(2, 3), function(x) quantile(x, 0.975))


# population proportions for each position
pi_est <- apply(rstan::extract(posterior)$pi,c(2,3),mean)
rownames(pi_est) <- issue_labels

# calculate importance by worse work situation for lower bound
ordinal_response_table <- matrix(NA,N_items,4)
for (j in 1:N_items) ordinal_response_table[j,] <- prop.table(xtabs(svyclean$weight~ordinal_responses[,j]))

# calculate total weight of importance for each issue, by level of worse work situation
importance_lower <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est_lower[j,] * (g-1)))
    for (p in 1:4) importance_lower[j,g] <- importance_lower[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est_lower[j,p]-psi_est_lower[j,])))
  }
}


# calculate importance by worse work situation for upper bound
ordinal_response_table <- matrix(NA,N_items,4)
for (j in 1:N_items) ordinal_response_table[j,] <- prop.table(xtabs(svyclean$weight~ordinal_responses[,j]))

# calculate total weight of importance for each issue, by level of worse work situation
importance_upper <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est_upper[j,] * (g-1)))
    for (p in 1:4) importance_upper[j,g] <- importance_upper[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est_upper[j,p]-psi_est_upper[j,])))
  }
}


#Mean importance
importance <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est[j,] * (g-1)))
    for (p in 1:4) importance[j,g] <- importance[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est[j,p]-psi_est[j,])))
  }
}

importance_table <- data.frame(issue=issue_labels,
                               importance_control=importance[,1],
                               importance_lower_control=importance_lower[,1],importance_upper_control=importance_upper[,1],
                               importance_treat=importance[,2],
                               importance_lower_treat=importance_lower[,2],importance_upper_treat=importance_upper[,2])


importance_table_print <- data.frame(
  Importance_Control = round2(importance_table$importance_control),
  Interval_control = paste0("(",round2(importance_table$importance_lower_control),"-",round2(importance_table$importance_upper_control),")"),
  Importance_Treat = round2(importance_table$importance_treat),
  Interval_Treat = paste0("(",round2(importance_table$importance_lower_treat),"-",round2(importance_table$importance_upper_treat),")")
)

rownames(importance_table_print) <- beta_treatment_teable$issue
colnames(importance_table_print) <- c("Not SVP", " ", "SVP", " ")

#app_imp_treat_betas.tex (in Latex file)
print(xtable(importance_table_print,
             caption = "Importance estimates given a worsened work situation or no change wit 95% posterior intervals."
             ,digits=3),
      type="latex",
      table.placement="!h", comment=FALSE, file = "Appendix_Table_A19.tex"
)

# ------------------------------------------------------------------------------
# Left-Right orientation: issue importance   -----------------------------------
# ------------------------------------------------------------------------------
#Issue Importance for Left-Right Orientation
load("Surveydata_fullsample.Rdata")
load("Posterior_fullsample_by-leftright.Rdata")

beta_chain <- rstan::extract(posterior)$beta 
N_sim <- dim(beta_chain)[1]
beta_est <- apply(beta_chain,c(2,3),mean)
beta_se <- apply(beta_chain,c(2,3),sd)
rownames(beta_est) <- issue_labels
colnames(beta_est) <- c( "Treatment (economic shock)")

# calculate which coefficients are significantly different from the mean coefficient for that variable
beta_mean_chain <- apply(beta_chain,c(1,3),mean)
beta_sign_chain <- array(NA,dim=dim(beta_chain))
for (sim in 1:N_sim) for (j in 1:N_items) beta_sign_chain[sim,j,] <- beta_chain[sim,j,] > beta_mean_chain[sim,]
beta_p <- apply(beta_sign_chain,c(2,3),mean)
beta_sig <- -1*(beta_p < 0.05) + (beta_p > 0.95)
beta_star <- array(c("(-)","","(+)")[beta_sig+2],dim(beta_sig))
beta_print <- data.frame(round(beta_est,3),beta_star)
beta_print <- beta_print[rev(inv_sort_order),]

#Confidence intervals beta coefficient
beta_est <- apply(beta_chain,c(2,3),mean)
beta_int <- t(apply(beta_chain,c(2),quantile,c(0.025,0.975)))

beta_treatment_teable <- data.frame(issue=issue_labels,
                                    beta_treatment=beta_est,
                                    beta_treatment_lb=beta_int[,1],beta_treatment_ub=beta_int[,2])

beta_treatment_table_print <- data.frame(
  Beta_Treatment = round2(beta_treatment_teable$beta_treatment),
  Interval = paste0("(",round2(beta_treatment_teable$beta_treatment_lb),"-",round2(beta_treatment_teable$beta_treatment_ub),")"),
  Significance = beta_star
)
rownames(beta_treatment_table_print) <- beta_treatment_teable$issue

# add average across issues
beta_treatment_table_print[N_items+1,c(2)] <- ""
beta_treatment_table_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_treatment_table_print)[N_items+1] <- "Average"
colnames(beta_treatment_table_print)[c(3)] <- ""

# add average across issues
beta_print[N_items+1,c(2)] <- ""
beta_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_print)[N_items+1] <- "Average"
colnames(beta_print)[c(2)] <- ""


# For interpretation: Calculate multiplicative effect from importance ----------
beta_multi <- -1*(1-exp(beta_est))
percent_matrix <- matrix(NA, nrow = nrow(beta_multi), ncol = ncol(beta_multi))

for (i in 1:nrow(beta_multi)) {
  for (j in 1:ncol(beta_multi)) {
    percent_matrix[i, j] <- paste0(signif(beta_multi[i, j] * 100, digits = 2), "%")
  }
}

multi_average <- paste0(round((-1 * (1 - exp(-0.085))) * 100, digits = 1), "%")
percent_matrix <- c(percent_matrix, multi_average)


# add the percent_matrix as a new column to the data frame
df <- cbind(beta_treatment_table_print, percent_matrix)
colnames(df) <- c( "Beta coefficients for treatment", "Interval", "", "Multiplicative effect")

#app_imp_treat_betas_multi.tex (in Latex file)
print(xtable(df,
             caption = "Coefficient estimates for variation in importance given a righ-wing political orientation. Coefficients
significantly higher than or lower than the average coefficient for the treatment variable across all issues are marked
with (+) or (-), respectively. The multiplicative effect is the natural exponent of the beta coefficient, indicating the impact of being treated on the disutility weight attached to voter-candidate disagreements."
             ,digits=2),
      type="latex",
      table.placement="!h", comment=FALSE, file = "Appendix_Table_A20.tex"
)

# Calculate Confidence intervals for beta 
confidence_level <- 0.95
beta_est_lower <- apply(beta_chain, c(2, 3), function(x) quantile(x, (1 - confidence_level) / 2))
beta_est_upper <- apply(beta_chain, c(2, 3), function(x) quantile(x, 1 - (1 - confidence_level) / 2))

## positions of alternatives
psi_est_lower <- apply(rstan::extract(posterior)$psi, c(2, 3), function(x) quantile(x, 0.025))
psi_est_upper <- apply(rstan::extract(posterior)$psi, c(2, 3), function(x) quantile(x, 0.975))


# population proportions for each position
pi_est <- apply(rstan::extract(posterior)$pi,c(2,3),mean)
rownames(pi_est) <- issue_labels

# calculate importance by worse work situation for lower bound
ordinal_response_table <- matrix(NA,N_items,4)
for (j in 1:N_items) ordinal_response_table[j,] <- prop.table(xtabs(svyclean$weight~ordinal_responses[,j]))

# calculate total weight of importance for each issue, by level of worse work situation
importance_lower <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est_lower[j,] * (g-1)))
    for (p in 1:4) importance_lower[j,g] <- importance_lower[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est_lower[j,p]-psi_est_lower[j,])))
  }
}


# calculate importance by worse work situation for upper bound
ordinal_response_table <- matrix(NA,N_items,4)
for (j in 1:N_items) ordinal_response_table[j,] <- prop.table(xtabs(svyclean$weight~ordinal_responses[,j]))

# calculate total weight of importance for each issue, by level of worse work situation
importance_upper <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est_upper[j,] * (g-1)))
    for (p in 1:4) importance_upper[j,g] <- importance_upper[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est_upper[j,p]-psi_est_upper[j,])))
  }
}


#Mean importance
importance <- matrix(0,N_items,2)
for (j in 1:N_items) {
  for (g in 1:2){
    issue_weight_tmp <- exp(sum(beta_est[j,] * (g-1)))
    for (p in 1:4) importance[j,g] <- importance[j,g] + ordinal_response_table[j,p]*sum(ordinal_response_table[j,]*abs(issue_weight_tmp*(psi_est[j,p]-psi_est[j,])))
  }
}


importance_table <- data.frame(issue=issue_labels,
                               importance_control=importance[,1],
                               importance_lower_control=importance_lower[,1],importance_upper_control=importance_upper[,1],
                               importance_treat=importance[,2],
                               importance_lower_treat=importance_lower[,2],importance_upper_treat=importance_upper[,2])


importance_table_print <- data.frame(
  Importance_Control = round2(importance_table$importance_control),
  Interval_control = paste0("(",round2(importance_table$importance_lower_control),"-",round2(importance_table$importance_upper_control),")"),
  Importance_Treat = round2(importance_table$importance_treat),
  Interval_Treat = paste0("(",round2(importance_table$importance_lower_treat),"-",round2(importance_table$importance_upper_treat),")")
)

rownames(importance_table_print) <- beta_treatment_teable$issue
colnames(importance_table_print) <- c("Pol. Left", "", "Pol. Right", "")

#app_imp_treat_betas.tex (in Latex file)
print(xtable(importance_table_print,
             caption = "Importance estimates given a worsened work situation or no change wit 95% posterior intervals."
             ,digits=3),
      type="latex",
      table.placement="!h", comment=FALSE, file = "Appendix_Table_A21.tex")
      
      
# ------------------------------------------------------------------------------
# Subgroup Analyses: High/low Environmental Attitudes and High/Low Income ------
#                    and issue importance as function of economic affectedness -
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# Above median environmental concern         -----------------------------------
# ------------------------------------------------------------------------------

rm(list= ls()[!(ls() %in% c('round2','issue_labels','issue_labels_signed', 'N_items'))])
load("Surveydata_envhigh.Rdata")
load("Posterior_envhigh_by-treatment.Rdata")

beta_chain <- rstan::extract(posterior)$beta 
N_sim <- dim(beta_chain)[1]
beta_est <- apply(beta_chain,c(2,3),mean)
beta_se <- apply(beta_chain,c(2,3),sd)
rownames(beta_est) <- issue_labels
colnames(beta_est) <- c( "High EA")

# calculate which coefficients are significantly different from the mean coefficient for that variable
beta_mean_chain <- apply(beta_chain,c(1,3),mean)
beta_sign_chain <- array(NA,dim=dim(beta_chain))
for (sim in 1:N_sim) for (j in 1:N_items) beta_sign_chain[sim,j,] <- beta_chain[sim,j,] > beta_mean_chain[sim,]
beta_p <- apply(beta_sign_chain,c(2,3),mean)
beta_sig <- -1*(beta_p < 0.1) + (beta_p > 0.90)
beta_star <- array(c("(-)","","(+)")[beta_sig+2],dim(beta_sig))
beta_print <- data.frame(round(beta_est,3),beta_star)

# add average across issues
beta_print[N_items+1,c(2)] <- ""
beta_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_print)[N_items+1] <- "Average"
colnames(beta_print)[c(2)] <- ""
colnames(beta_print)[c(1)] <- "High EA"

# collapse beta chain into four categories
helper1 <- apply(beta_chain[,1:4,],1,mean) #mean value per row (=respondent) for issue 1-4
helper2 <- apply(beta_chain[,5:8,],1,mean)
helper3 <- apply(beta_chain[,9:12,],1,mean)
helper4 <- apply(beta_chain[,13:16,],1,mean)

cov1 <- c(helper1, helper2, helper3, helper4)
covariate1 <- matrix(cov1,nrow=3000,ncol=4,byrow=FALSE)
betcol_chain <- sapply(list(covariate1), I, simplify="array")

betcol_est <- apply(betcol_chain,c(2,3),mean)
betcol_se <- apply(betcol_chain,c(2,3),sd)
rownames(betcol_est) <- c("Environment", "Environment/Economy", "Economy", "Other")
colnames(betcol_est) <- c("High EA)")


betcol_mean_chain <- apply(betcol_chain,c(1,3),mean)

betcol_sign_chain <- array(NA,dim=dim(betcol_chain))
for (sim in 1:N_sim) for (j in 1:4) betcol_sign_chain[sim,j,] <- betcol_chain[sim,j,] > betcol_mean_chain[sim,]
betcol_p <- apply(betcol_sign_chain,c(2,3),mean)
betcol_sig <- -1*(betcol_p < 0.05) + (betcol_p > 0.95)
betcol_star <- array(c("(-)","","(+)")[betcol_sig+2],dim(betcol_sig))
betcol_print <- data.frame(round(betcol_est,3),betcol_star)

# add average across issues
betcol_print[4+1,c(2)] <- ""
betcol_print[4+1,c(1)] <- round(apply(betcol_mean_chain,2,mean),3)
rownames(betcol_print)[4+1] <- "Average"
colnames(betcol_print)[c(2)] <- ""

betcol_print1 <- betcol_print

# ------------------------------------------------------------------------------
# Below median environmental concern         -----------------------------------
# ------------------------------------------------------------------------------
rm(list= ls()[!(ls() %in% c('betcol_print1', 'round2','issue_labels','issue_labels_signed', 'N_items'))])
load("Surveydata_envlow.Rdata")
load("Posterior_envlow_by-treatment.Rdata")

beta_chain <- rstan::extract(posterior)$beta 
N_sim <- dim(beta_chain)[1]
beta_est <- apply(beta_chain,c(2,3),mean)
beta_se <- apply(beta_chain,c(2,3),sd)
rownames(beta_est) <- issue_labels
colnames(beta_est) <- c( "Low EA")

# calculate which coefficients are significantly different from the mean coefficient for that variable
beta_mean_chain <- apply(beta_chain,c(1,3),mean)
beta_sign_chain <- array(NA,dim=dim(beta_chain))
for (sim in 1:N_sim) for (j in 1:N_items) beta_sign_chain[sim,j,] <- beta_chain[sim,j,] > beta_mean_chain[sim,]
beta_p <- apply(beta_sign_chain,c(2,3),mean)
beta_sig <- -1*(beta_p < 0.1) + (beta_p > 0.90)
beta_star <- array(c("(-)","","(+)")[beta_sig+2],dim(beta_sig))
beta_print <- data.frame(round(beta_est,3),beta_star)


# add average across issues
beta_print[N_items+1,c(2)] <- ""
beta_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_print)[N_items+1] <- "Average"
colnames(beta_print)[c(2)] <- ""
colnames(beta_print)[c(1)] <- "Low EA"

# collapse beta chain into four categories
helper1 <- apply(beta_chain[,1:4,],1,mean) #mean value per row (=respondent) for issue 1-4
helper2 <- apply(beta_chain[,5:8,],1,mean)
helper3 <- apply(beta_chain[,9:12,],1,mean)
helper4 <- apply(beta_chain[,13:16,],1,mean)

cov1 <- c(helper1, helper2, helper3, helper4)
covariate1 <- matrix(cov1,nrow=3000,ncol=4,byrow=FALSE)
betcol_chain <- sapply(list(covariate1), I, simplify="array")

betcol_est <- apply(betcol_chain,c(2,3),mean)
betcol_se <- apply(betcol_chain,c(2,3),sd)
rownames(betcol_est) <- c("Environment", "Environment/Economy", "Economy", "Other")
colnames(betcol_est) <- c("Low EA")


betcol_mean_chain <- apply(betcol_chain,c(1,3),mean)

betcol_sign_chain <- array(NA,dim=dim(betcol_chain))
for (sim in 1:N_sim) for (j in 1:4) betcol_sign_chain[sim,j,] <- betcol_chain[sim,j,] > betcol_mean_chain[sim,]
betcol_p <- apply(betcol_sign_chain,c(2,3),mean)
betcol_sig <- -1*(betcol_p < 0.05) + (betcol_p > 0.95)
betcol_star <- array(c("(-)","","(+)")[betcol_sig+2],dim(betcol_sig))
betcol_print <- data.frame(round(betcol_est,3),betcol_star)

# add average across issues
betcol_print[4+1,c(2)] <- ""
betcol_print[4+1,c(1)] <- round(apply(betcol_mean_chain,2,mean),3)
rownames(betcol_print)[4+1] <- "Average"
colnames(betcol_print)[c(2)] <- ""

betcol_print2 <- betcol_print


# ------------------------------------------------------------------------------
# Above median income                        -----------------------------------
# ------------------------------------------------------------------------------
rm(list= ls()[!(ls() %in% c('betcol_print1', 'betcol_print2', 'round2','issue_labels','issue_labels_signed', 'N_items'))])
load("Surveydata_inchigh.Rdata")
load("Posterior_inchigh_by-treatment.Rdata")

beta_chain <- rstan::extract(posterior)$beta 
N_sim <- dim(beta_chain)[1]
beta_est <- apply(beta_chain,c(2,3),mean)
beta_se <- apply(beta_chain,c(2,3),sd)
rownames(beta_est) <- issue_labels
colnames(beta_est) <- c( "High income")

# calculate which coefficients are significantly different from the mean coefficient for that variable
beta_mean_chain <- apply(beta_chain,c(1,3),mean)
beta_sign_chain <- array(NA,dim=dim(beta_chain))
for (sim in 1:N_sim) for (j in 1:N_items) beta_sign_chain[sim,j,] <- beta_chain[sim,j,] > beta_mean_chain[sim,]
beta_p <- apply(beta_sign_chain,c(2,3),mean)
beta_sig <- -1*(beta_p < 0.1) + (beta_p > 0.90)
beta_star <- array(c("(-)","","(+)")[beta_sig+2],dim(beta_sig))
beta_print <- data.frame(round(beta_est,3),beta_star)


# add average across issues
beta_print[N_items+1,c(2)] <- ""
beta_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_print)[N_items+1] <- "Average"
colnames(beta_print)[c(2)] <- ""
colnames(beta_print)[c(1)] <- "High income"

# collapse beta chain into four categories
helper1 <- apply(beta_chain[,1:4,],1,mean) #mean value per row (=respondent) for issue 1-4
helper2 <- apply(beta_chain[,5:8,],1,mean)
helper3 <- apply(beta_chain[,9:12,],1,mean)
helper4 <- apply(beta_chain[,13:16,],1,mean)

cov1 <- c(helper1, helper2, helper3, helper4)
covariate1 <- matrix(cov1,nrow=3000,ncol=4,byrow=FALSE)
betcol_chain <- sapply(list(covariate1), I, simplify="array")

betcol_est <- apply(betcol_chain,c(2,3),mean)
betcol_se <- apply(betcol_chain,c(2,3),sd)
rownames(betcol_est) <- c("Environment", "Environment/Economy", "Economy", "Other")
colnames(betcol_est) <- c("High income")


betcol_mean_chain <- apply(betcol_chain,c(1,3),mean)

betcol_sign_chain <- array(NA,dim=dim(betcol_chain))
for (sim in 1:N_sim) for (j in 1:4) betcol_sign_chain[sim,j,] <- betcol_chain[sim,j,] > betcol_mean_chain[sim,]
betcol_p <- apply(betcol_sign_chain,c(2,3),mean)
betcol_sig <- -1*(betcol_p < 0.05) + (betcol_p > 0.95)
betcol_star <- array(c("(-)","","(+)")[betcol_sig+2],dim(betcol_sig))
betcol_print <- data.frame(round(betcol_est,3),betcol_star)

# add average across issues
betcol_print[4+1,c(2)] <- ""
betcol_print[4+1,c(1)] <- round(apply(betcol_mean_chain,2,mean),3)
rownames(betcol_print)[4+1] <- "Average"
colnames(betcol_print)[c(2)] <- ""

betcol_print3 <- betcol_print

# ------------------------------------------------------------------------------
# Below median income                        -----------------------------------
# ------------------------------------------------------------------------------
rm(list= ls()[!(ls() %in% c('betcol_print1', 'betcol_print2','betcol_print3','round2','issue_labels','issue_labels_signed', 'N_items'))])
load("Surveydata_inclow.Rdata")
load("Posterior_inclow_by-treatment.Rdata")

beta_chain <- rstan::extract(posterior)$beta 
N_sim <- dim(beta_chain)[1]
beta_est <- apply(beta_chain,c(2,3),mean)
beta_se <- apply(beta_chain,c(2,3),sd)
rownames(beta_est) <- issue_labels
colnames(beta_est) <- c( "Low income")

# calculate which coefficients are significantly different from the mean coefficient for that variable
beta_mean_chain <- apply(beta_chain,c(1,3),mean)
beta_sign_chain <- array(NA,dim=dim(beta_chain))
for (sim in 1:N_sim) for (j in 1:N_items) beta_sign_chain[sim,j,] <- beta_chain[sim,j,] > beta_mean_chain[sim,]
beta_p <- apply(beta_sign_chain,c(2,3),mean)
beta_sig <- -1*(beta_p < 0.1) + (beta_p > 0.90)
beta_star <- array(c("(-)","","(+)")[beta_sig+2],dim(beta_sig))
beta_print <- data.frame(round(beta_est,3),beta_star)


# add average across issues
beta_print[N_items+1,c(2)] <- ""
beta_print[N_items+1,c(1)] <- round(apply(beta_mean_chain,2,mean),3)
rownames(beta_print)[N_items+1] <- "Average"
colnames(beta_print)[c(2)] <- ""
colnames(beta_print)[c(1)] <- "Low income"

# collapse beta chain into four categories
helper1 <- apply(beta_chain[,1:4,],1,mean) #mean value per row (=respondent) for issue 1-4
helper2 <- apply(beta_chain[,5:8,],1,mean)
helper3 <- apply(beta_chain[,9:12,],1,mean)
helper4 <- apply(beta_chain[,13:16,],1,mean)

cov1 <- c(helper1, helper2, helper3, helper4)
covariate1 <- matrix(cov1,nrow=3000,ncol=4,byrow=FALSE)
betcol_chain <- sapply(list(covariate1), I, simplify="array")

betcol_est <- apply(betcol_chain,c(2,3),mean)
betcol_se <- apply(betcol_chain,c(2,3),sd)
rownames(betcol_est) <- c("Environment", "Environment/Economy", "Economy", "Other")
colnames(betcol_est) <- c("Low income")


betcol_mean_chain <- apply(betcol_chain,c(1,3),mean)

betcol_sign_chain <- array(NA,dim=dim(betcol_chain))
for (sim in 1:N_sim) for (j in 1:4) betcol_sign_chain[sim,j,] <- betcol_chain[sim,j,] > betcol_mean_chain[sim,]
betcol_p <- apply(betcol_sign_chain,c(2,3),mean)
betcol_sig <- -1*(betcol_p < 0.05) + (betcol_p > 0.95)
betcol_star <- array(c("(-)","","(+)")[betcol_sig+2],dim(betcol_sig))
betcol_print <- data.frame(round(betcol_est,3),betcol_star)

# add average across issues
betcol_print[4+1,c(2)] <- ""
betcol_print[4+1,c(1)] <- round(apply(betcol_mean_chain,2,mean),3)
rownames(betcol_print)[4+1] <- "Average"
colnames(betcol_print)[c(2)] <- ""

betcol_print4 <- betcol_print

# ------------------------------------------------------------------------------
# Appendix Table A22: Combine subgroup results into table          -------------
# ------------------------------------------------------------------------------

# Merge the first variable of each data frame
merged_df <- cbind(betcol_print1$High.EA., betcol_print2$Low.EA, betcol_print3$High.income, betcol_print4$Low.income)

# Rename the columns of the merged data frame
colnames(merged_df) <- c("High EA", "Low EA", "High income", "Low income")
rownames(merged_df) <- c("Environment", "Environment/Economy", "Economy", "Other", "Average")


print(xtable(merged_df,
             caption = "Coefficient estimates for variation in importance given a worsened work situation in different subgroups: high and low environmental attitudes (EA), high and low income."
             ,digits=3),
      type="latex",
      table.placement="!h", comment=FALSE, file = "Appendix_Table_A22.tex", include.rownames = FALSE
)


#End of File